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Fourier Series: Review 


Fourier Series: Review 


A function or signal a(t) is called periodic with period T if x(t + T) = x(t). All “typical” 
periodic function x(t) with period T can be developed as follows 
Equation: 


ao a 27 : 27 
ys S k 
x (t) 5 + a a;,COs ( F t) + b;sin (: T 7 


where the coefficients are computed as follows: 


Equation: 
T/2 
a-zf (t}o0s (it) a Ce, Leone’) 
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The natural interpretation of [link] is as a decomposition of the signal x(t) into individual 
oscillations where a, indicates the amplitude of the even oscillation cos (k22t) of frequency k/T 
(meaning its period is T’//k), and by indicates the amplitude of the odd oscillation sin (k-2 t) of 
frequency k/T. For an audio signal x(t), frequency corresponds to how high a sound is and 


amplitude to how loud it is. The oscillations appearing in the Fourier decomposition are often also 
called harmonics (first, second, third harmonic etc). 


Equation: 


Note: one can also integrate over [0, T] or any other interval of length T’. Note also, that the 
average value of x(t) over one period is equal to ag/2. 


Complex representation of Fourier series 
Often it is more practical to work with complex numbers in the area of Fourier analysis. Using the 
famous formula 
Equation: 
e!* = cos (a) + jsin (a) 
it is possible to simplify several formulas at the price of working with complex numbers. Towards 


this end we write 
Equation: 


1,. fx 1 — 4 
acos (a) + bsin (a) = a= (e +e) + bas (e® —eJ*)= 3 (a — bjje* + 3 (a + bj)je 7 


i) 


From this we observe that we may replace the cos and sin harmonics by a pair of exponential 
harmonics with opposite frequencies and with complex amplitudes which are conjugate complex to 
each other. In fact, we arrive at the more simple complex Fourier series: 

Equation: 


. ; Lf | 
z(t) = S> XpeF2th/T with X,=— / ip (t)e—s2mkt/T gy 
k=—00 T 0 


Note that X is complex, but x is real-valued (the imaginary parts of all the terms in }> X,¢e27*4/T 


add up to zero; in other words, they cancel each other out). The absolute value of X; gives the 
amplitude of the complex harmonic with frequency k/T (meaning its period is T’/k); the argument 
of X;, provides the phase difference between the complex harmonics. If x is even, X; is real for all 
k and all harmonics are in phase. 


To verify [link] note that by [link] and [link] we have for positive k 
Equation: 


x= = [ve (t)[cos(2rkt/T) — jsin(Qrkt/T)|dt = 5 (a — bj). 


For negative k we note that X_, = X ‘ by [link], where ()° denotes the conjugate complex. By 
[link], the X; are exactly as they are supposed to be. 


Properties 
¢ Linearity: The Fourier coefficients of the signal z(t) = ca(t) + y(t) are simply 


Equation: 


Zr = cX,+ Y, 


¢ Change of frequency: The signal z(t) = x(At) has the period T/A and has the same Fourier 
coefficients as x(t) — but they correspond to different frequencies f: 


Equation: 


Cc [o-@) 
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T/A 
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The equation on the right allows to read off the Fourier coefficients and to establish 2; = X x. (For 
an alternative computation see Comment 1) 


Comment 1 
Equation: 


1 T/r Xr T/r 
Z, = a/ z (the?) dt = z/ ataje re at 
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1 n —2njks/T 
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¢ Shift: The Fourier coefficients of z(t) = x(t + d) are simply 
Equation: 
q.= Xperts 
The modulation is much more simple in complex writing then it would be with real coefficients. 
For the special shift by half a period, i.e., d = T'/2 we have Z, = X,e!™* = (1X, 


¢ Derivative: The Fourier series of the derivative of x(t) with development [link] can be 
obtained simply by taking the derivative of [link] term by term: 


Equation: 


oe 2m) Cg 
x! (t) = 3 Xp be + ePehlt 


k=—oo 


Short: when taking the derivative of a signal, the complex Fourier coefficients get multiplied by 
hams. Consequently, the coefficients of the derivative decay slower. 


Examples 


e The pure oscillation (containing only one real but two complex frequencies) 
Equation: 


z(t) = sin (2zt/T) X,=-X_1= i 


or B, = 1/2 = —B_,, or b; = 1, and all other coefficients are zero. This formula can be 
obtained without computing integrals by noting that 
sin (a) = (e* — e~J) / (27) = (7/2) (e 9% — e%) and setting a = 2nt/T. 

e Functions which are time-limited, i.e., defined on a finite interval can be periodically 
extended. Example with T’ = 27: 
Equation: 


tl for @<t<7 4 
z(t)= -1 for —-7<t<0 by = 2Be = — for odd k>1 
c fort=0,7 e 
and all other coefficients zero. Note that c is any constant; the value of c does not affect the 


coefficients b,. We have for0 <t <7 
Equation: 


4 1 1 
1 = —{ sin (t) + —sin (3t) + =sin (5t)+.. ) 
T 3 5 


Note that for t = 0 the value of the series on the right is 0, which is equal to x(0—) + 2(0+), 
the middle of the jump of x(t) at 0, no matter what c is. Similar for t = 7. 


Discrete Fourier Transform 


Discrete Fourier Transform 


The Discrete Fourier Transform, from now on DFT, of a finite length sequence 
(x0, ..-, £K_1) is defined as 
Equation: 


MM. = ne %*F (k=0,..,K—1) 


To motivate this transform think of z,, as equally spaced samples of a T’-periodic 
signal x(t) over a period, e.g., 2, = x (nT/K). Then, using the Riemann Sum 
as an approximation of an integral, i.e., 


Equation: 

K-1 T 

nT T 
f—->: =-2 f (t)dt 

~~. kK Kk 0 
we find 
Equation: 

K-1 T 


0 


Note that the approximation is better, the larger the sample size K is. 


Remark on why the factor K in [link]: recall that X; is an average while x; is a 
sum. Take for instance k = 0: Xo is the average of the signal while xo is the sum 
of the samples. 


From the above we may hope that a development similar to the Fourier series 
[link] should also exist in the discrete case. To this end, we note first that the DFT 
is a linear transform and can be represented by a matrix multiplication (the 
“exponent” 7’ means transpose): 

Equation: 


T ‘a 
LO, --> LK] = DFT x : (x0, ep O41) : 


The matrix DF'T’x possesses K lines and K rows; the entry in line k row n is 
e 2nskn/K Aw few examples are 
Equation: 


DFT, =(1) DFT, = eg 


j -l-3 


S 

y 

a 

| 
ee eS 


The rows are orthogonal[footnote] to each other. Also, all rows have 
length[footnote]/K . Finally, the matrices are symmetric (exchanging lines for 
rows does not change the matrix). So, the multiplying DFT with its conjugate 
complex matrix (DF Tx) we get K times the unit matrix (diagonal matrix with 
all diagonal elements equal to K’). 

The scalar product for complex vectors x = (#1, Z9,..., 2) and 


y = (Y1, Y2;---, yk) is computed as 
Equation: 


L Y= L1Yy + LQYQ +.» TLKYR; 


where a’ is the conjugate complex of a. Orthogonal means z - y = 0. 
Length is computed as 


2 2 
2) =Waa = aj0, + to, +..+2Kt_% = x1 +|x2 +...4|ex|? 


Inverse DFT 
From all this we conclude that the inverse matrix of DFT’; is 


IDFT x = (1/K)-(DFTx) . Since (e~*) = e® we find 
Equation: 


rper rr (n = 0,...,K —1) 


Spectral interpretation, symmetries, periodicity 


Combining [link] and [link] we may now interpret x; as the coefficient of the 
complex harmonic with frequency k/T in a decomposition of the discrete signal 
Xn; its absolute value provides the amplitude of the harmonic and its argument 
the phase difference. 


If x is even, x, is real for all & and all harmonics are in phase. 


Using the periodicity of e?" we obtain x, = 2, when evaluating [link] for 
arbitrary n. Short, we can consider x, as equally-spaced samples of the T- 
periodic signal x(t) over any interval of length T: 

Equation: 


K/2-1 
t, = L pe 
n=—K/2 


—Injkz 


Similarly, x, is periodic: x, = 2.4K. Thus, it makes sense to evaluate x, for any 
k. For instance, we can rewrite [link] as 
Equation: 


1 KPA 


L,=— Le 
K n=—K/2 


Qnjkz 


Since x; corresponds to the frequency k/T’, the period K of x, corresponds to a 
period of K’/T in actual frequency. This is exactly the sampling frequency (or 
sampling rate) of the original signal (K samples per T' time units). Compare to 
the spectral repetitions. 


However, the period T of the original signal z is nowhere present in the formulas 


of the DFT (cpre. [link] and [link]). Thus, if nothing is known about 7;, it is 
assumed that the sampling rate is 1 (1 sample per time unit), meaning that kK = T’ 


FFT 


The Fast Fourier Transform (FFT) is a clever algorithm which implements the 
DFT in only K log (£) operations. Note that the matrix multiplication would 
require K? operations. 


Matlab implements the FFT with the command fft( x) where z is the input 
vector. Note that in Matlab the indices start always with 1! This means that the 
first entry of the Matlab vector z, i.e. z(1) is the sample point x9 = x (0). 
Similar, the last entry of the Matlab vector z is, i.e. x(4’) is the sample point 
Qe =2((k —-1)T/K)=a(T-—T/kK). 


Fourier Integral 


Fourier Integral 


Continuous-time signals which are not periodic can still be understood 
as superpositions of pure oscillations where now all frequencies are 
present in the signal. The coefficients of the oscillations can be 
computed as follows: 

Equation: 


The representation as a superposition takes then the following form: 
Equation: 


We call the Fourier transform of and write also instead of 

to indicate clearly which signal has been transformed. The “Fourier 
spectrum”, or simply the spectrum, or also the “power spectrum” of the 
signal is the squared amplitude . This is the function usually plotted, 
while the phase of _ is not shown. Nevertheless, the plots are usually —and 
erroneously— labeled with instead of (see [link]). 


A signal is called bandlimited if its Fourier transform is zero for high 
frequencies, i.e. for large . Similarly we say that a signal is time-limited if 
it is zero for large times, i.e., for large . By Heisenberg's principle a 
bandlimited signal can not be time-limited. Since bandlimited signals are of 
great importance, there is a need to study signals which are not time-limited 
and, thus, the Fourier integral. 


Properties 


e Linearity: 
Equation: 


e Convolution 
Equation: 


e Change of time scale 
Equation: 


e Translation in time and frequency 
Equation: 


e Symmetries and Fourier pairs The symmetry of [link] and [link] leads 
one to consider and as a Fourier pair. Indeed, the Fourier 
transform of isalmost : . Clearly, the 
symmetry is not perfect since is in general complex, while is 
real. However: If is symmetric, i.e. then is 
real-valued, and vice versa! 


In summary: Symmetric real signals have symmetric real Fourier transforms 
and vice versa. As we will see below, they also possess the same energy. 


Energy and Power 


Energy and Power 


The energy of a continuous-time signal x(t) is given as 
Equation: 


lei? = f n? (t)dt 


Plancherel's theorem says (for more information see Comment 2): If the signal 
z(t) has finite energy then its Fourier transform X(f) has the same energy: 
Equation: 


Ico = fixe ag= fo 2? @ae= Hel? rieiemerey eas 


Comment 2 Plancherel theorem is a result in harmonic analysis, first proved by 
Michel Plancherel. In its simplest form it states that if a function f is in both 

Ly (IR) and Lz (IR), then its Fourier transform is in Dy (IR); moreover the 
Fourier transform map is isometric. This implies that the Fourier transform map 
restricted to Ly (IR) Lz (IR) has a unique extension to a linear isometric map 
D2 (IR) > Lz (IR). This isometry is actually a unitary map. 


Periodic signals have of course infinite energy; therefore, one introduces the 
power P,, of the signal x(t), which is the average energy over one period. Energy 
is measured in Joule, power is measured in Watt=Joule/Sec. 


The analog of Plancherel's theorem is Parseval's theorem which applies to T- 
periodic signals and says 
Equation: 


1 T/2 ee) 
P= PS / x (t)|?dt = X,,|" periodic case 
ii a (t)| 2 | | 


* 


We may derive Parseval's theorem as follows, using [link] and Jal? =a*a: 
Equation: 
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— z[ eora= 2 S > XpePHIT at 
7 T J_7/2 T J_r/2 " 


k=—0o 
ok 


i or S- Qakt/T S- Qnnt /T 

— — Xe? Xe? oe dt 
T —-T/2  k=—00 n=—0o 
I hacks S- 2nkt/T 3 * 5 —jannt /T 

= — X,e7"" Ase dt 
T a a ee nN=—0o 


00 0 «] T/2 
2 Sy ae / eiem(k—n)t/T ay 
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2S * 
= > XX}, = > IX," 


k=—00 k=—0o 


ej2n(k—n)t/T 


Here, we used that equals T’ when k = n (since e® = 1), but 


equals 0 when k # n since (since e/* =cos (s) + 7 sin (s), which are integrated 
over several periods). 


A similar computation can be carried out for Plancherel's equation. However, 
some difficulties arise due to the integrals over infinite intervals (see Comment 3 
below). Also, a justification of Plancherel could be given by performing a limit of 
infinite period in Parseval's equation (see Comment 4 below). 


For finite discrete signals the analog is simply the fact, that DFT is unitary up to a 
stretching factor. More precisely, the matrix DFT’ leaves angles intact and 
stretches length by VK. Intuitively, one may think of the DFT as a rotation and a 
stretching. In other words, to perform a DFT simply means to change the 
coordinate system into a new one, and to change length measurement by a factor 


VK . Thus: 


Equation: 


= a t,°=— iscrete case 
KS K? : Kk * 


Note that the DFT Fourier coefficients are complex numbers; thus, the absolute 
value has to be taken (for a complex number a we have |a| * =a-a’, whichis 
usually different from a?—unless a is by chance real valued). 


Comment 3 A “hand-waving” argument for Plancherel's theorem runs as follows, 
using [link] and |a|? = a-a’: 


Equation: 
[emda [f° Po xinerrar a 
= / - / _X (Hem hap / _ X (ge? "dg ‘at 
/ - / _X (eas / : X(g) e "dg dt 
FF [somes [oma 
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/ _X(HX(A'at -/ - IX (f)/ae = [|X]? 


Co Jamt(f—9) ty — = 
Thereby, the step —. e€ dt = 6(f — g) would require some more care, 
but we content ourselves here with this intuitive computation. 


An important example of a Sinc and the un- The power 
band-limited signal: the normalized version in spectrum of 
sinc-function. comparison. sinc is Rect(f). 


Rect(f) 


sin(ax)(xx) 


\eewwa ewews beewe feewd SWewe GOewe COeee Geewe Seee ewS POETS EET | 1 1 L 1 1 1 1 
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Comment 4 With Parseval's equation established one may provide a derivation of 
Plancherel's formula which is more convincing, tho more difficult to make 
rigorous. Assume that the signal x(t) is time-limited, say defined on 

—T/2<t < T/2. Its Fourier-series coefficients are 


7 T/2 
Kpeiyt se , 


thus T - X; = X (k/T). We may interpret x as being periodically extended and 
use Parseval's equation which says P,» = _|_X4|*; this allows the following 
computation 

Equation: 


a (t)e 27*/T dt. Comparing to the Fourier transform we find 


T/2 0 | 
I|a|[? = / ln (t)|"dt =T-P,=T >> X, = S| ST? |X;! 
fe) 


=e k=—0o ee | 
[© @) 1 (ee) 

= > Fxanr~ [ ix@Pae 
k=—0o — oo 


In the last step we use a Riemann sum approximation of the integral. It remains to 
extend this to signals which are not time-limited. The approximations made are 
more believable than the step using the Dirac delta (see Comment 3 above). 
Although more intuitive and believable, these arguments are harder to make 
rigorous than making the computation of Comment 3 rigorous by passing through 
a correct, direct computation of the integrals avoiding the Dirac delta. 


Examples 


Examples 


e The pure oscillation (containing only one frequency) 
Equation: 


sin (2nt/T) {sin (2nt/T)} (f) = 2(6(-1) -6 (0) 


This formula can be obtained without computing integrals by noting that 

sin (x) = (e%” — e~*) / (23) = (j/2) (e-* — e%”). Its power is P = 1/2 (see [link]). 
The perfect low (frequency) pass function: 

Equation: 


1 if —1/2<f<1/2 
0 else. 


sinc (t) := F {sinc(t)} (f) = Rect (f) = { 


and more general (to pass exactly the frequencies f €|—f., fe) 
Equation: 


Qfe-sinc(2f.t) F{2fpsinc (2f-t)} (f) = Rect ( if ) _ {* if —f.<f< fe 


2fe 0 else. 
Both, x and X are symmetric and real. Plancherel's formula allows to compute the energy of the sinc: 
Equation: 
fore) 1/2 
sinc I? =f x?(pag= frag =1. 
—oo —1/2 


More generally, the energy of 2f, - sinc (2f-¢) amounts to 2f,. Note that the sinc is not time-limited; it can't 
be by the Heisenberg principle since it is bandlimited. 

¢ The Dirac function 6(¢) is often symbolically written as 
Equation: 


oo” if t=0 


5 = {0 oe FHM =A 


Clearly, the Dirac function is not really a function, and it has infinite energy. However, most manipulations 
work fine also for 6(¢). This is again an illustration of the Heisenberg principle. The Dirac function is the 
extreme case which is sharply located in time, but has no characteristic frequency (all frequencies are present 
with equal strength). The properties of the Dirac function are best understood in terms of integrals: 


Equation: 
. 0) if 0 
[owswae= (oO 5a ba 


0 else 


As a special case, the convolution with a function g is again g: 
Equation: 


co [o-e} 


S(a—t)g()dt= f 5(t)g(a-t)dt = g(a) 


—oo 


{5%9} (a) = | 


—oo 


short: 6*g(t) = g(t). For the shifted Dirac function 6, (t) = 6 (t — b) we have 
Equation: 


CO CO 


5(t—Bg(a— sar = [ Solo --b)idsSale—9) 


—oo 


{60} (a) = | 5.()g(a- eat = [ 
—oo —oo 

short, convolution with 6, produces a shift by b: 6,*g (t) = g(t — b). 

Double exponential: 

Equation: 


2 


z(ty=e "4 X(f)= 14 4qafe 


Note that x and its Fourier transform X are real and symmetric. The power spectrum is 
|X (f)|? =4/ (1+ 41? f?) * Since a is not differentiable at 0, the Fourier transform X decays somewhat 


slowly: high frequencies are quite strong in this signal in order to make the sharp peak at 0. With this 
example, we may compute the energy directly: 


Equation: 
co 2 Co 
le“? = / (e#) dt = 2 | "tae "|" =i, 
—0o 0 


One-sided Exponential: 
Equation: 


e' ift>0 1 
t = xX = 
au) {f else. (f) 1+ 2njf 


The Fourier transform X is complex with power spectrum |.X (f)|? = 1/ (1 + 4m? f?). Since a is not even 
continuous at 0, the Fourier transform X decays even slower than for the double exponential: high 
frequencies are even stronger in this signal in order to make the jump at 0. With this example, we may verify 
Plancherel's theorem: 


Equation: 
co 2 co 
ai? = f (e~*) at = | Mat =~ ale = 1/2. 
0 0 
Equation: 
1 , aa 1 1 1 
cag = fo egg = ge arctan On|, = 5 (n/2= (=m /2)) = 1) 


The Gaussian Kernel is practically invariant under the Fourier transform (see Comment 5) 
Equation: 


e(t)=e"? — X(f) = VIme Cr"? 


Here, it is easy to verify ||x|| = ||X|| via a substitution t = 27f. The computation is somewhat harder and 
yields ||a||? = 7 ~ 1.7725. 


Comment 5 From Probability theory, we know that (see “characteristic function of a Gaussian distribution”) 
Equation: 


[etetnas _ ef /2 
V 20 


Now replace f by 27f and multiply with V 27 to find the Fourier transform. For the energy: 


Equation: 


oo 7 2 = oo 
lel? =f (e°?)ae= vee feta = va 


TT lo) 


1 
Vr 
variance 1/2, and thus integrates to 1. 
The Mexican Hat (also called Ricker Wavelet in Geophysics) is the negative second derivative of the 
Gaussian: 
Equation: 


t 


where we use in the last step, that —>e~ * constitutes the probability density of a Gaussian variable with 


a(t) = (1 = Pee X(f) = An? /2nf eo Onf)?/2 


Both, the Gaussian kernel and the Mexican hat are very useful since they are well located both in space and in 
frequency (see [link]), meaning that the main portion of their energy stems from a narrow range of locations 
as well as a narrow range of frequencies. Thus, the may be used as low-pass, respectively band-pass filters. 
The Mexican hat is a wavelet; wavelets are used to determine which frequencies contribute the main portion 
of the energy at a specific time. To this end, a wavelet needs to be well localized in time as well as in 
frequency. 

The Dirac Comp (peigne de Dirac) of step 7 

Equation: 


tA, (t):=7 >) 6(t—nr) F{r-A,}(f) = 3 6(f —k/r) 


k=—00 


To verify this formula, we choose the integration interval [—7/2, 7/2] and write 
Equation: 
1 7/2 os 7/2 a. 7/2 sh 
XE = =| r+ A, (t)e P™/ "dt = By 6(t —nr)e P™/ "dt = / 5 (t)e Pm!" dt = 1 
= n —7/2 


T 7/2 —7/2 


The Fourier-representation becomes 
Equation: 


r-A, (t) _ SS edenkt/r 


k=—00 


The Gaussian 
kernel. Its Fourier 
transform has the 

identical shape. 


The Mexican hat is 
the negative second 
derivative of the 
Gaussian kernel. 


Fourier transform 
of the Mexican hat. 


This formula will be crucial in the reconstruction formula in the Nyquist-Shannon Theorem. Note: One 


should not try to evaluate 7 - A, (t) ; indeed, it is not really a function, since it is formed by Dirac terms. 
Even its power is infinite. 


Estimation of Spectrum and Power via DFT 


Estimation of Spectrum and Power via DFT 
Spectrum for finite energy signals 


The Fourier integral is intended for non-periodic signals, since the periodicity will lead to problems 
with the convergence of the integral. In fact, the Fourier integral [link] is defined rigorously only for 
signals with finite energy (see next section below). Recall that finite energy signals can never be 
periodic. For an alternative way of defining X(f) for periodic signals see the next bullet just below). 


However, in order to make the DFT also useful for non-periodic finite energy signals, we can consider 
any signal x(t) as being restricted to the interval |—L/2, L/2] and then extended periodically. Denote 
this new L-periodic signal by % (t) with Fourier coefficients X yz. When L is large, comparing [link] 
with [link] gives 

Equation: 


2s 1 pe ; i. - 7? , 1 k 
Ea z(t —jamkt/L yy pee / t —jamkt/L yy a5 ae 
: Z| p20 rg cs LE 


—oCo 


The notation X x instead of X, should remind the reader, that x itself is not periodic. 


If not the entire signal z but only some K samples z, = x (nL/K) of the signal are given, we may 
use the DFT as in [link] and combine it with [link] to get fork = —K/2,..,K/2-—1 


Equation: 
a k > os 
finite energy: x = L-Xpm 2p - 


Note that this approximation improves in quality as L becomes larger; the factor L/K indicates that 
would have to be made proportionally larger at the same time. For an illustration see [link], [link], and 
[link], left column, with L = 20, K = 50. Imagine that we reduce L to half, i.e. to L* = 10, leading 
to a new signal x. This would mean to discard half of the samples, thus leading to K” = 25. Since the 


discarded samples are practically zero, using [link] shows that LE Lox. This agrees with [link] since 
L* /K* = L/K and X (k/L”) = X (2k/L). 


signal with 50 sample points | DFT | / (Sampling rate) 
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Spectral Estimation via DFT. A finite energy signal (Gaussian 
kernel). 


signal with 30 sample points | DFT | / (Sample size) 
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Spectral Estimation via DFT. The sampled periodic signal 
z(t) = 1+ sin (27t)+ sin (3 - 27t)+ sin (4 - 2nt) (period T = 1, 
length L = 3). The larger the sample-size A, the better are the 
locations (0, 1, 3 and 4) and amplitudes (1, 0.5, 0.5 and 0.5) of the 
delta peaks of X approximated by |%;|/K. 


signal with 50 sample points 


‘ F i 
1 ots if ll ¢ 
a Ii ee 
Wale NIV ALY 
oly) sft 
wi OVP OV 
-i} | | ; 
-1 0 1 


| DFT | / (Sample size) 
1 if 


0.5 


Le 


Spectral Estimation via DFT. The sampled periodic signal 
z(t) = 1+ sin (27t)+ sin (3 - 27t)+ sin (4 - 2nt) (period T = 1, 
length LZ = 3). The larger the sample-size A, the better are the 
locations (0, 1, 3 and 4) and amplitudes (1, 0.5, 0.5 and 0.5) of the 
delta peaks of X approximated by |%;|/K. 


4o 0 0 


Spectrum for periodic signals 


For T-periodic signals, the Fourier integral [link] does not converge in the usual sense. One can still 
define a Fourier transform X(f) consisting of Dirac delta functions: 
Equation: 


X(f) = 3) Xa5 (f —k/T) ~ > aed (F —/T), 
k k 


With this setting we have x(t) = f X (f)e??"/*df also for periodic signals. Note that the Delta- 
functions in [link] are “infinitely large” at f = k/T, in agreement with [link] where L should ideally 
be “infinitely large”. 


However, “infinitely large” values are not useful when one wants to read off an amplitude. Comparing 
[link] and [link] we see that the DFT x; should be normalized differently for non-periodic finite energy 
signals than for periodic signals. Dividing the DFT by the sampling rate for non-periodic signals 
provides an estimation of the Fourier transform via [link]. Dividing the DFT by the sample size for 
periodic signals provides an estimation of the peak heights of the Dirac delta functions in [link] (see 
[link] and [link)). 


Energy and Power from Sampling 


First of all, note that the energy of a periodic signal is infinite because there are infinitely many periods 
with the same energy. On the other hand, the power of a finite energy signal is zero since it is the 
average energy of an infinitely long interval. 


Short, only one measure is meaningful for a signal: either the power or the energy, depending on 
context. Nevertheless, we may compare power and energy of the various ways of representing a signal: 
as a finite energy signal, as a signal on a given finite interval (extending it periodically) or through its 
samples. 


Comparing the power of a 7’-periodic signal with its samples over an interval of length LD, where L is 
an integer multiple of the period T, we find in the frequency domain (using [link] and [link]) 
Equation: 


ee) K/2 K/2 1 ; 
P; (cont) = > eel? = S> bel = S- K2 |z;.| = ae (sampled) 
k=—0o k=—K /2 k=—-K/2 


Computing in the time domain (using a Riemann Sum) we arrive at the same conclusion: 
Equation: 


1 L/2 ; 1 K/2 L L 2 1 K/2 ; 
(cont) = + / 1p Ola = = > ale ng 2 DF | acicied) 


Comparing the energy for finite energy signals and for their samples on an interval of length L we find 
Equation: 


Le 
Energy = ||z||? ~ / |x (t)|’dt = L- P, (cont) ~ L- P, (sampled) 
2 


where we used [link] in the last step. 


In summary, using sufficiently many samples taken over sufficiently large intervals should allow to 
compute the power or energy of signals in either of the three representations, thereby approximating 
power with energy per time over a long interval. 


As an important conclusion we note that the power of a periodic signal should not depend on the 
number of samples taken, at least approximatively, and as long as sufficient many samples are taken. 
We will later make precise how many samples are sufficient. Similar for the energy of a finite energy 
signal. 


Anti 
As a final remark we note that the name “power spectrum” used for |X (f)|”, | Xx|" and |Zx.|° is 
appropriate since they indicate how the power is distributed of the spectrum of frequencies, i.e., which 
frequencies contribute much and which less to the power. 


Sampling: Review 


Sampling: Review 


Let us consider a signal x(t) which is uniformly sampled meaning that the 
samples are x, = x(nT) (n=... — 2,—1,0,1, 2, ...). The sampling frequency 
f. = 1/r is the number of samples per time unit. 


In practical applications, the signal will be sampled only over an interval, say of 
length Z and into K samples. Then 7 = L/K as for discrete signals. 


Spectral copies (also called spectral repetitions or images) 


The sampled signal can be written as the Dirac comb of step 7 modulated by 


x(t) 


Equation: 


The factor 7 is added to obtain a Fourier transform that does not depend on T 
and to preserve an averaging property (the integral of x, provides a good 
approximation of the integral of x since fx, (t)dt = )>°°._., ra (nr) is the 


Riemann sum approximating the fz (t) dt as asum of rectangles of width 7.) 


Using [link] we may write 
Equation: 


[o.@) 


te (t)=2(t)>7A,() = S° a (t) - ef27h/r 


k=—oo 


Note, that [link] is not a Fourier representation of x, (t) (the “coefficients” x(t) 
of the sinusoids depend on ¢ instead of k). Computing the Fourier transform of 
the sampled signal x, (t) once more, using [link], we find 

Equation: 


CO CO 


Xf) = > {eW@e® rb (p= So x(F- 4/7) 


k=—0o k=—0o 


We observe that the Fourier transform, and also the spectrum of the sampled 
signal is periodic with period f. = 1/7, the sample rate (see [link] right). 
However, we see now that X, is composed of copies of X, shifted by a 
multiple of the period, i.e., shifted by k/7 = k - f.. These are called spectral 
copies. 


Sampling and the discrete Fourier transform 


Finite energy signals: Using [link] again we find the Fourier transform of the 
sampled signal x, (t) 
Equation: 


[e,6) 


Xe(f) = > ra(nr) {d(t-nr)}(f) =7 9 Ly eam int 


n=— CO n=— CO 


Setting now f = k/L and recalling r = L/K we get approximatively 
Equation: 


K/2 
finite energy: X.(k/L) ~ + S- Bye PUKLIML/K) — 7 Ry, 
n=—K/2 


This agrees with [link] for finite energy signals and is actually valid for all k, 
since both sides are periodic with period K. However, it is only useful if x is of 
finite energy, similar to [Link]. 


For periodic signals we could write X and X, as sums of Dirac delta functions 
as in [link]. Doing so, [link] confirms that also periodic signals possess spectral 
copies when sampled. 


For periodic functions, the relation [link] shows that the terms 7x, will attempt 
to approximate the Dirac-shape of X,, meaning that the values of 7; will be 
huge for k/L close to a peak location but very small otherwise (see [link], 


[link], and [link], compare discussion after [link]). In order to identify the 
amplitude X;, of the Dirac delta functions better, one uses the earlier 
approximation 

Equation: 


periodic: Ap & —S]7 eg. 


Sampling at 300 
dpi appears to be 
sufficient to keep 
the quality of the 
image. 


Sampling at 50 dpi 
introduces aliasing. 


Sampling of a 
detail at S500dpi 
reveals high- 
frequency content. 
These frequencies 
will leak from the 
spectral copies into 
the low frequencies 
of the main spectral 
copy when sampled 
at too lowa 
frequency. 


Aliasing 


Assume that the continuous-time signal x(t) is bandlimited, meaning there is 
some B > 0 such that X(f) = 0 for |f| > B. If the sample rate had been too 
small, namely 

Equation: 


fe 2B 


then this copies would overlap and the original signal can not be recovered. The 
sampled signal shows artifacts called aliasing (recouvrement); these artifacts 
manifest is erroneous low frequency content. Such content spills or leaks from 
the spectral copies into the main spectral period (see [link]). 


Nyquist Shannon sampling theorem: 


Assume that the continuous-time signal x(t) is bandlimited, meaning there is 
some B > O such that X(f) = {a(t)}(f) =0 whenever |f| > B. 


If this signal x is sampled uniformly at a sufficient rate, namely at f. > 2B, 
then x(t) can be mathematically reconstructed perfectly from only those 
discrete samples. 


Reconstruction (Ideal low pass filtering) 


If the sample rate has been high enough though, namely 
Equation: 


Jo 28 
then, the Fourier transform of the original signal can be recovered by cutting off 


the periodic copies: 
Equation: 


where f, has to be chosen such that B < f. < f, — B (see [link] right). (The 
first copy to be cut off lies over the interval | f. — B, f. + B].) A possible 
choice is f. = f,/2. 


In the time domain, this corresponds to convolving the sampled signal with the 
ideal low-pass filter (sinc). More precisely: Translating [link] into the time 
domain, then using the linearity of the convolution and [link], we find the ideal 
low-pass filtering (reconstruction) formula: 

Equation: 


x(t) Le (t)*2f_sinc (2f.t) 


= € y x (nr)d(t — nr)*2f-sinc (2 ft) 
=: (27,) » x (nr)sinc (2f, (t — nr)) 


Note that the pre-factor 27 f, = 2f-/f- takes into account the mismatch 
between sampling rate f, and the cut-off frequency f,. When choosing 


f. = fe/2, it disappears (becomes 1). 


Practical reconstruction: Recall that the sinc filter is not realizable since it 
requires using all samples of the past and the future. Since a bandlimited signal 
can not be time-limited at the same time (Heisenberg), this would in theory lead 
to infinitely many operations. 


If fe >> 2B, then a filter different from the ideal sinc-filter can be used. In 
fact, the reconstruction x (t) = x, (t)*b(¢) is valid for any low-pass filter b 
(not only for b(t) = 2f,sinc (2f,t)) as long as we assure that the spectrum of 
the filter b is equal to 1 for all frequencies in |—B, B] and zero for frequencies 
larger f. — B. (Note that we do not specify the phase). Usually, one designs a 
filter such that 


= 


. The filter's spectrum equals 1 in an interval |—f, fp| called passing band. 
It must contain the signal's frequency-band |—B, B). 

2. The filter attenuates all frequencies beyond f, — fo. 

3. The transition band of the filter is [f,, fe — fp]. This range of frequencies 

can be chosen as large as |B, fe — B] by setting f, = B and as short as 

desired, but always containing the point f,/2 by setting f, smaller but as 

close to f,/2 as desired. 


The Dirac The spectrum of the sampled signal (note the spectral 
Comb copies) and a practical reconstruction filter: It's passing 

modulated band contains the band |—B, B] of the signal and it 
bya attenuates all spectral copies. 


signal. 


The low-pass FIR filter of 
length 61 and ideal cutoff at 
fe=1/8 =0, 125 as 
produced by Matlab. The 
filter values are in red. For 
comparison in blue the 
corresponding samples of 
the ideal filter 
2f-sinc (2f,t) at 


t = —30,..,—1,0,1,..30. 


digital filter in black (sample points in red); sinc for comparison (blue) 


transition 
band 


; assing band ; ' attenuation 
1 P § e——¥ 
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band 


DFT of the signals on the left. The actual 
transition band of the filter is [0. 1, 0. 15] 
. Note that the sinc-samples provide a 
filter with sharper cutoff but slightly 
lower quality in the pass-band since it 
requires infinitely many samples. The 
power computed from either side (recall 
[link]) is 
P, = 0.003891 = 0. 237'4/61 ~ 2f,/ 
length. 


IDFTI? of filter (black and red); IDFTI? of sinc for comparison (blue) 


All this assures that we keep the entire spectral information of x(t) contained in 
[—B, B], yet remove all spectral copies (see [link]). Sampling an ideal 
reconstruction filter is a method to a design a low pass filter. Here, point 2 in 
the above list helps avoiding aliasing in such a filter. However, the resulting 
filters are not ideal as FIR filters possible since not all samples can be used, 
resulting in Gibbs-like phenomena (see [link] on the right). 


Matlab uses the command fir to compute finite impulse response filters 
(FIR). A larger transition band allows for shorter filters to be used, thus 
speeding up computation. [link] shows a length 61 FIR filter with ideal cutoff 
frequency at f, = 1/8 = 0.125. Its transition band is roughly [0. 10, 0. 15]. As 
is common in digital filtering, the sampling frequency is assumed to be f, = 1. 


Power and Energy from Sampling 
We summarize a few facts. 


Good quality when above Nyquist We are now in a position to make more 
precise how many samples are enough to well approximate the power of a 
signal via [link], or its energy per time over a large interval via [link], 
depending on the context. Roughly, the signal must be sampled at least at 
Nyquist rate so that the samples and their DFT faithfully represent the signal 
and its Fourier transform. Clearly, sampling at a much higher rate than Nyquist 
will improve the approximation. 


Power independent of sampling rate We conclude that changing the sampling 
rate will not change power, as long as we stay above Nyquist, and at least 
approximately. Before we study how to change the sampling rate in the next 
sections, we give a quick demonstration of [link] using a simple band-limited 
signal: a filter. 


Ideal filter and FIR filter The energy of the ideal filter with cutoff frequency f, 
is ||2f-sinc (2f-t)||? = 2+ f, which follows easily from the power spectrum 
being 1 on an interval of length 2 - f, and zero else, using Plancherel. To 
compute the power of a digital low pass FIR filterb = (bo, ..., 6x1) of length 
K with approximate cutoff frequency f, we study its DFT b: we note that 

b, = 1 for roughly 2, of the indices k and b; = 0 for the other indices. 
Note that f. = 1 and0 < f, < 1/2 for a digital filter. Using [link] we find 
Equation: 


2 1 
= qa tie = 2f-/K = 2fc/L. 


This fits perfectly with [link] since for a digital filter f. = 1 and, thus, K = L. 
The approximation becomes better, the sharper the transition band of the filter b 
, L.e., the longer the filter is. 


Compare again with [link] for a concrete filter with length K = 61, 
fe = 1/8 = 0.125 and power P; = 0. 003'891. 


Decimation and Downsampling 
Decimation and Downsampling 


Pure Downsampling 


The easiest case consists in reducing the sampling rate by simply dropping samples. This procedure is 
called Downsampling. However, we need to be careful on the effect of this procedure. We note 
immediately, that the new sampling rate fg = f./M needs to be still above Nyquist, ie., f. > M-2B 
in order to avoid aliasing. 


¢ Downsampling Assume that f. > M - 2B. 
Equation: 


(Yo. Y1sY2---) > {M > (yo, ym; Yom,---) _ (Zo, 21, Z2,---) 


Effect in the frequency domain In continuous time, downsampling the signal y(t) corresponds to 
passing to the signal z(t) = y(Mt). Indeed, sampling z at the same frequency as y we obtain the 
samples 2, = yxy as we Should. In summary, using the properties of the Fourier transform 
Equation: 


z(t)=y(Mt) Z(f)= a (az) 


This tells us, that the spectral copies of Z, (the Fourier transform of the sampled signal z, (¢)) should 
have the same overall shape and at the same distance f, as those of Y,, but they are stretched in 
frequency by a factor M and squeezed in amplitude by a factor M. 


Indeed, we arrive at the same conclusion computing via the sampled signals, using [link]: 
Equation: 


; 1 1 
Z. (f) = Voraer™ _ S> ryace Panter te S> (Mr) ype 27 F/ MAM) _ — (4) 
; ; M ; M M 


Recall that Yq (f), the Fourier transform of y sampled at frequency fy = f-/M, looks like Y, except 
that its spectral copies lie at distance fg, thus M times closer than those of Y.. 


For the discrete Fourier Transform we should expect to see roughly the same behavior. Recall, though, 
that the relation between continuous and discrete Fourier transform is only approximative. 


downsampled digital filter in black (sample points in red); sinc for comparison (blue) IDFTI? of downsampled filter (black and red); IDFTF of sine for comparison (blue) 
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The result of downsampling the filter from [link] by a factor M = 2 
on the left, and its power spectrum |DFT|? on the right (black and 
red). Again, the ideal sinc filter in blue for comparison. Note that the 
filter's cutoff frequency has increased by M (stretching of the 
spectrum by M) and that its |DFT|” values have decreased by 
1/M? = 0. 25 (see text). The power computed from either side is 
0.003829 (recall [link]): downsampling of a signal with 2B < f./M 
does not change power. 


Aliasing: We give two examples. 


First, downsampling to a sampling rate that is too low can lead to aliasing. Downsampling the image of 
[link] left leads to the same effect as visible in the same figure center. 


Second, downsampling a simple signal such as the filter b of [link] by M = 2 results in a new filter b’ 
with twice the cutoff frequency, i.e. f{ = 1/4 = 0. 25, but with a value b, = 1/2 = 1/M over the 
pass-band, thus a power spectral value 1/4 = 1/M? over the pass-band (see [link]). No aliasing occurs 
since the condition f, > M - 2B is satisfied. 


Power: We consider first the case of a T-periodic signal y. Then, z(t) = y(Mt) has period T/M. 
Substituting s = tM with ds = Mdt we get 
Equation: 


1 T/M ; mM tM ; 1 rt : 
3S: t)|“dt = — Mt)|"dt = — ds = P 
auf Ora=s Pf wampa= sf w(oPas=e, 


From the simple properties we know Z;,| a Yel , (same complex amplitudes[footnote], but 
T/M —- 
Ee 
belonging to different frequencies). Using Parseval we see again that the power does not change. See 
properties of Fourier Series. 


There is no contradiction between Z;, = Y;, and Z(f) = 47Y (+). Both express that z(t) = y(Mt), 


both allow to conclude that power does not change under downsampling, and both imply that the DFT 
of zis M times smaller and corresponds to frequencies which are M times further apart than those of y. 
Indeed, using [link] we get 

Equation: 


K 1 1. 
— £, = — KY, = —tk 


7h) f=afr ~ Mf M M 


ft 
and using [link] with D equal to the correct period and number of samples (T'/M and K/M for z(t) 


respectively T and K for y(t)) we get again 
Equation: 


; K/M k K 1 1 k 1 K(k 1. 
2h pote eee | = ae! Va ee | ae ee) 
7 T/M°\T/M) TM \MT/M) MT \T) M 


fat 


Note that the energy of a finite energy signal would change under downsampling: computing in time 
Equation: 


la? =f lewptae= fo wenn itar= Ff Iy(s)Pas = ill 


A simliar computation can be done in frequency. See Comment 6. 


Comment 6 Computing in frequency with f = Mg anddf = Mdg: 
Equation: 


Z|? = [za = [|r (2) 


For a discrete signal we are naturally in the periodic case. From the above we should expect that power 
stays approximatively the same under downsampling provided that f, > M - 2B. Also, we noted 
earlier that power should not depend on the sampling rate, as long as the samples faithfully represent the 
signal, and at least approximatively. 


a=, f Wolds = FIrIP 


For the simple signal 6, the filter from [link], we may verify this explicitly. Denote the downsampled 
filter by b’ (see [link] for an illustration with M = 2). Since no aliasing occurs during downsampling, 


the pass-band is now M times longer, meaning that M times more of the b x are different from zero 
a (2 
b! " are by M2 smaller 


(this makes the power decrease by M7“). Finally, the sample length is now M times shorter (this 
increases the power my MV; recall [link]). 


(this makes the power increase by M). Further, their power spectral values 


All in all, power is not changed, at least approximatively. One finds Py = 0. 003’829 which has to be 
compared to the power P; = 0. 003’891 of the original filter. 


To obtain a low-pass filter one would have to normalize b’ to b” = M - b’. 


Decimation 


Let us now drop the assumption f, > M - 2B. 


To resample a signal x(t) at an M times lower rate, a first attempt would be to discard all but every M/- 
th sample: z, = py = x (kM). This step is called downsampling as we have seen above. However, 
to avoid aliasing effects caused by downsampling below Nyquist rate a low-pass filtering at cutoff 

B* = f./ (2M) is required before downsamling. The filter used is called anti-aliasing filter. The new 
utilized bandwidth will be only B’. 


The procedure of first applying an anti-aliasing filter and then downsampling is called decimation. 


Agreeably, the filtering before downsampling destroys information about x(t). However, this loss 
occurs in a controllable manner: It removes high-frequency information. For an audio signal, we loose 
the high pitch sound. For an image we loose sharpness of edges. This should be compared to an 
uncontrolled loss of quality when no anti-aliasing filter is applied (see [link]). 


In summary, Decimation by M means to resample at M times lower rate fg = f-/M and consists of 
two steps: 


1. low-pass filtering the samples at cutoff frequency sir 


Equation: 


(xo, 21,2 )> nn : + ( ) 
Or Ula W2y--s 2M Yo, Y1,Y2;--- 


2. Down-sampling (in French: decimation) 
Equation: 


(Yo, Yt, Y2 . =) = 1 M- (Yo, YM> Y2M) - . ai) — (20; 21,22; . -) 
Low-pass filtering will reduce power, as high frequencies are attenuated. Down-sampling will leave the 
new power roughly the same. 


The matlab commands are decimate and downsamp1le. 


Interpolation and Upsampling 


Interpolation and Upsampling Interpolation 


Let us now look at increasing the sample rate. Since it is less obvious how to achieve this, let us first 
consult theory. 


Given are the samples z,, x nt of aband-limited signal x t taken at frequency f. T is 
above the Nyquist rate JB. We want to compute the samples z; of z t at a higher sampling rate 
fu N  f,-. This means that the new sampling step shouldbe ff, feN 7 IN Ot, 

zr avkrn. 


Since the original sampling rate f. B is above Nyquist, we can in theory reconstruct the entire 
signal  t using the reconstruction formula [link] using f. jf. .Once the continuous-time (finite 
energy) signal z t is obtained, we only need to sample it att kr Nk f,. This reads as follows 
Equation: 


Ze Zz kr x kv Ln +. Ln 


This formula allows indeed to compute zz from zy, at least in principle. However, a closer look at 
theory is required to understand the effect when using only finite many samples. The reconstruction 
formula [link] is best understood in the frequency domain: it amounts to removing the spectral copies 
of X. f via filtering with cut-off frequency f. fe .To make this filtering step visible we need to 
write [link] in form of a convolution. To this end, we write 

Equation: 


2 
=| 3 


m 
2k Ym —— Ym W N Um 


where the new sequence y,,, is obtained by “upsampling” and is given as: 
Equation: 


Ba mN n 


om m N 


The convolution [link] allows for more convenient data processing via digital filtering and for a 
simple spectral interpretation. 


Interpolation by N or resampling at N times larger rate consists of the following steps: 


1. Up-sampling (in French: “interpolation”) 
Equation: 


N N N 


2. Multiplication by N and Low-pass filtering at cutoff frequency —y using the ideal filter 


— ahs 
N N° 
Equation: 
m 
yyy IN xa Ym W a ae 
digital filter upsampled by N=3 in black; sinc/N for comparison (blue) IDFTI? of filter upsampled by N=3; IDFTF of sinc/N for comparison (blue) 
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The result of upsampling the filter from [link] by a factor NV on 
the left, and its power spectrum on the right (black and red). 
Again, the ideal sinc filter (divided by N to adjust for the zero- 
samples) in blue for comparison. Note that the filter's cutoff 
frequency has decreased by N (contraction of the spectrum by NV) 
and that its values have not changed (see text). The power 
computed from either side is 0.001'297 and has decreased by a factor 
N from the original power of 0.003'829: only upsampling of a signal 
with Bf. decreases power by the upsampling factor. 


Spectral picture of Interpolation by NV 


In analogy to the decimation we set z t xz t N .Then, the samples of z ¢ taken at the same rate 
f- constitute the samples of x taken at rate f,,. Analogously to the decimation we find quickly 
Equation: 


zt «tN Zf NXfN 


which indicates how the spectrum at rate f, Nf, is obtained: the spectrum at rate f, gets 
contracted in the frequency axis by N and expanded in amplitude by NV. 


Note that the spectral copies of Z, are at distance f, just like those of X¢. 


interpolated digital filter in black (sample points in red); sinc for comparison (blue) uIDETE of interpolated filter (black and red); IDFT? of sinc for comparison (blue) 
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The result of interpolating (upsampling and filtering) the filter from 
[link] by a factor NV on the left, and its power spectrum 
on the right (black and red). Again, the ideal sinc filter (no need to 
divide by N since the zero-samples are now corrected through the 
filtering step) in blue for comparison. Note that the fundamental 
period contains now only one spectral copy, as it should, as a result 
from the filtering. The spectrum is contracted by N in frequency and 
expanded in amplitude by NV . The power computed from either side 
is 0.003'891 and almost identical to the original power of 0.003'829: 
interpolation of asignal with Bf, does not change power. 


We may break the procedure down into the individual steps: 


(i) Upsampling (introducing the zero-samples) leaves the Fourier transform, and thus the spectrum 
almost intact, leading only to a rescaling of the frequencies (contraction of X): Indeed, the Fourier 
transform Y, of the samples y,,, becomes 

Equation: 


ee i TYme j ufmr ran,e ITNT X, Nf 


m n 


The corresponding signal y t (with samples y,, at sampling rate f,) is not of interest. If you want to 
know about it any way, see Comment 7. Note, however, that the fundamental period of Y. amounts to 
f- and contains N copies of X Nf (see [link]). 


Comment 7 For clarity: the Fourier transform of y is found by removing the spectral copies of Y. 
outside f. fe. These copies are caused by sampling. The Fourier transform of y consists of NV 


contracted copies of X at distance f, N of each other which are caused by upsampling. Using the 
reconstruction formula [link] with the samples y,,,, sample rate f, and correct pass-band f. yields 
the signal y t Ym funn . Using that while m for all integer 
m we find quickly that the samples of yt areindeedy kr = yy, i.e., 

x £ a 


(ii) Multiplication with N restores the average value of the samples. The Fourier transform is now 
NX. Nf which consists of copies of NX fN ,orZ f , at distance f. N (as for Y., there are N 
copies in one period). 


(iii) The digital low-pass filtering of Ny, at cut-off frequency N removes all of the copies of 
NX fN except the ones centered at 0, fe, fe etc. and leaves only one copy per period, in other 
words, only copies at distance f.. What remains is exactly Z. f as we have pointed out earlier. (see 
(link]). 


Power and Interpolation 


Similarly as with the decimation the power of a periodic signal does not change under interpolation. In 
fact, we only needtouse Z f © NX fN andreplace M by N inthe computation done with the 
decimation. We conclude that the power of discrete samples does not change under interpolation 
provided that f. B, at least approximatively. 


To verify this, let us move through the 3 steps above. Step (i), upsampling, reduces power by a factor 
N since the sum of squares of the samples is the same (the zeros added don't contribute), but there are 
now NN times more samples. Power is an average. 


Step (ii) obviously multiplies power with N . Step (iii), the low-pass filtering, removes VV 
spectral copies and leaves only 1, thus divides power by NV. All steps together leave the power as it is. 
Interpolation 


Interpolation illustrated in the Spectral Domain with NV 
Horizontal arrows indicate spectral copies. To keep the sample rate at the same value f, one 
changes from X,, resp. Yq to Z, (see d). 
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N 
inserting N 
zeroes contracts the 
Fourier transform 
by a factor N 


digital low-pass 
filtering with cut- 
off frequency at 
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Decimation 


Decimation illustrated in the Spectral Domain with M 
Horizontal arrows indicate spectral copies. To keep the sample rate at the same value f, one 
changes from Yq to Z, (see d). 
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Sampling Rate Conversion 


Sampling rate conversion 


In general, sampling rates are converted only in rational fractions such as 
from 32 to 48 kHz by a factor of 3/2 which are performed by appropriate 
sequences of up- and down-sampling. 


Practical considerations: 


e First interpolate, then decimate (preserve signal quality) 

e the low-pass filters used after upsampling (2nd step of interpolation) 
and before downsampling (1st step of decimation) can be combined 
into the one filter with the smaller pass-band. 

e Interpolation and Decimation by large factors should be done in steps, 
keeping NV and J/ below or at 6 at most. 

e Matlab commands: decimate, downsample, interp, 
upsample, resample. See the corresponding help manual (help 
decimate etc). 

e The ideal sinc filter is not realizable since it requires knowledge of the 
infinite past and future; one uses FIR filters with linear phase to avoid 
audible artifacts. 


Models of Noise 


Models of Noise 


White Noise In a nutshell, a sequence (€1, €2,€3,...) is called white noise if its entries €, are independent 
identically distributed random numbers. 


e “Independence” means that there is no information or relation between the members €;; in other words, 
knowing or observing the sequence until time ¢, i.e., (€1,...,€) does not allow to predict the next 
member €;41 any better than if nothing had been observed. 

e “Random” means that in principle, an entry can take any value in some given set with certain chance, and 
that it is are not in advance which value it will take. 

e “Tdentically distributed” means that each entry has the same chances to assume a possible value. 


White noise is an example of an stationary, ergodic series. Stationarity means that the statistics don't change 
over time. Ergodicity means that statistical measures such as mean and variance can be estimated by observing 
enough samples of one single sequence (€1, €2, €3,... ); the mean or expectation[footnote] of an entry E[e;| 
can be estimated as the sample mean (1/NV) (€1+...+e€w). 

Recall that IE[X] denotes the expectation of the random variable X. 


An example of a stationary, non ergodic series is the one where €; equals 1 or —1 with equal probability, and 
all other €% are equal to €;. Clearly, the sample mean of one single sequence is then either 1 or —1, but not 0 as 
it should. 


An example of white noise is the error €; introduced by quantization, i.e., €, = yr — ©, where x; isa 
“typical” signal before and y;, the signal after quantization. Check of Randomness and Independence: as we 
are observing the first t samples (€1,...,€,) we have no indication whatsoever on the quantization error of the 
t + 1st sample — unless the signal is very special. [An example of an atypical signal would be one that is 
already quantized: after observing 500 times an error 0 we start to suspect that the future errors will also be 0.] 
Check of the distribution: a quantization done by rounding to the third decimal, e.g., will result in errors that lie 
between [—0. 00049999. .. ,0. 0005] where all values in this interval are equally likely. This means, e.g., that 
€x is negative with chance 1/2 and that, e.g., €; is within [0. 0002, 0. 0003] with chance 1/10. Since this is the 
same of all entries €,, they are identically distributed. 


Spectral analysis of stationary signals and series 


By their own nature, similarly to periodic signals, stationary signals and series have no finite energy. As a 
simple example consider the sequence €; = +1 with random sign. Since ee = 1 for all k, the energy of the 
sequence is clearly infinite. In order to arrive at a meaningful spectral analysis one defines the power of a 
stationary signal x(t) as the time average of the energy: 

Equation: 
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Note that for a periodic signal definition [link] gives the same value as [link]. 


While periodic signals possess a natural Fourier expansion into a series, we need to take time-averaged, 
windowed Fourier transforms for stationary signals; also, due to randomness, one needs to take averages over 
different realizations in order to obtain a deterministic non-random spectral descriptor. Most useful is the 
power spectrumS( f) which corresponds to the square of the absolute value of the time-averaged windowed 
Fourier transform: 

Equation: 
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Sometimes it is (erroneously) written as S (f) = |E(f)|?. erroneously because it should be averaged: 
S(f) = E|E(f)|’. The famous Wiener-Khinchine theorem says, that 
Equation: 
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where r(k) = EE [én - €n—~] = IE [Ex - €0] is the auto-correlation of the process. Note S(f) = S(f +1). 


Connection to Discrete Fourier transform (DFT) In practice, a finite length signal €9,...,€%_1 is 
interpreted as periodically repeated (with period /’). Its Fourier transform is then a periodic series again, called 
DFT, which can be computed via an algorithm called FFT in only K log (/’) computations according to the 
formula 

Equation: 


K-1 
DFT: é, = S ape 
k=0 


Recall that Matlab starts indices always with 1, thus interpreting the first entry of the vector é as corresponding 
to frequency 0 (see matlab help fft). 

If we set K = 2D + 1 we find that €,, appears on the right side of [link]. If we assume that L is large so that 
we can neglect the limit, we find that S(f) indicates what to expect of the squared FFT coefficients (i.e., the 


FFT-power-spectrum) on average: 
Equation: 


Ble.) ~K-s() 


forn = 1,.., K. Vice versa, the samples of the power spectrum S(f) for f = k/K can be estimated from 


averaging the FFT coefficients over several noise signals each of the form €1,... , €ar+1: 
Equation: 
k 1 2 
S(—)2+ SE/é, 
x) =~ Ele! 
fork = —L,...,0,...,Z (it is more natural to represent S over a symmetric interval; recall that 
E_L = EL41,-++,€0 = E2L41). 


Flat spectrum of zero mean white noise 


The power spectrum $(f) of white noise with zero mean (IE[e;,| = 0) can be computed as follows. First we 
find 
Equation: 


r(0)=Elej] =o? and ry=EE[ex]-Eleo] =0 k#0 


by independence. By [link] we find 
Equation: 


S(f) =o? = Power 


for all f. Note that the power spectrum of white noise is flat. Its constant value is equal to the variance of the 
zero mean noise. The constant value must clearly be equal to the power. (Also: 
P, =1/K Se ~ E [e2] = o? from statistics.) 


Consequently, since the mean of the DFT, i.e. E\é, | is constant, the average of the DFT when computed over 
many realizations of the noise is nearly flat[footnote]; however, the FFT of any single realization will show 
large oscillations (see [link] 3rd). 

matlab demo: plot(mean(abs(fft(0.1*(rand(1000,1000)-0.5),[ ],2)..2)));axis([0 1000 0 0.01]) 


Nevertheless, the power can be estimated from one realization: using [link], using ergodicity to estimate the 
expected value IE]. . . | simply as the average over all samples, and using [link] we get 
Equation: 


This computation also confirms that the spectrum of white noise is flat and that its constant value is equal to 
the power. 


In fact, a direct computation from [link] shows that the FFT output €,, is a white noise as well, however with 
variance Ko? and so [link] (the first approximation in [link]) holds for white noise exactly, not only 
approximatively. 

Intuitively, the formula S (f) = o* means that all frequencies are present in white noise with equal overall- 
amplitude. The spectrum being flat is a direct consequence of the independence between the noise terms. We 
should recall, though, that strictly speaking the spectrum is flat only as an average over many noise 
realizations. 


Application: Quantization In many application, such as inference in wireless communication and such as in 
quantization in signal processing, a commonly accepted model of the effect of a source of errors on the signal 
z(t) is the so-called additive white noise model: 

Equation: 


Lp > +Ek > Yk 


For a quantization with precision A the error is the quantification error €, = y;, — £4 and it is determined by 
the signal itself. [One could argue, that €; is not random since it is completely determined once x; is known. 
Still, the model is useful since we can usually not predict €,41 from observing 71,.., £z.] 


For quantization using rounding (matlab: round|footnote]) it can be shown that the noise power per sample 


2 , : ee . 
amounts to Pouant.noise = a = . When using one more bit for quantization, then the error A is half as 


large, thus the power 4 times smaller, which amounts to roughly -6 dB. In other words, the power of the 
quantization noise is proportional to -6dB times the number of bits used. 


The error €, is in this case uniformly distributed on the interval [—A/2, 4/2]. 


Using the analog of Parseval's equation and recalling that S(f) is the power spectrum (analog of the square of 
the Fourier transform) we have indeed: 
Equation: 


A? 1 fe/2 A? 
Se (f) = SpReet(*) Pouant noise = fe [ py Se (f)df = 12. 


For K samples of noise taken over a time interval of length K’/f. we have, thus, approximatively 
Equation: 

1 K-1 1 K-1 A2 
= FP Oiisntncise = 2 


Note that the FFT increases power by K. The relation [link] becomes exact when taking expected values FE. 
The approximation improves the larger K is, since the left side is an estimator of the variance o”, which is 
A? /12 for quantization with precision A. 


Application: Interference A further example of a situation where an additive white noise proves useful is 
wireless transmission of a binary signal under interference. Here, bits may flip from 0 to 1 and vice versa since 
detection is not perfect. The error €, = yr — £» is here determined by the interfering signal, and the 
configuration of the decoder. Note that the possible values of each ¢; is either 0 (no flip) or 1 (flip). Without 
specific information on the interference, the chance of a flip is independent of the time k, and independent of 
the past occurrence of flips. Thus, white noise is a very reasonable model for €z. Clearly, the probability 

Ple,, = 0] will be close to 1 if only little inference is present and will decrease the stronger the inference. 


Gaussian noise As an important special case we mention Gaussian white noise, where the common 
distribution of the ¢, is Gaussian, or “normal”: This model assumption is standard whenever nothing is known 
on the distribution. It makes sense, e.g., as model for an overall error which is composed of several small 
unknown errors. (compare Central Limit Theorem) 


Colored Noise A sequence (1, 2, n3, ... ) is called colored noise if its terms n;, are random and possess a 
relation or dependence between them. Consequently, the power spectrum of colored noise is not flat, but 
possesses certain prevalent frequencies — hence the name “colored” (recall that the frequencies of light waves 
correspond to colors). 


One of the most simple ways to produce colored noise is to filter white noise. For instance, mz = Ex + €x41 is 
colored since the entries m, are no longer independent: mp = €9 + €; and m, = €; + €2 contain the same 
number € as an additive term. Similarly, nz, = €% — €%41 is colored. 


Adopting a continuous-time notation (for convenience) we write 
Equation: 


m(t) = e(t) + e(¢+ 1) n(t) = e(t) — e(t +1) 


with Fourier transform 
Equation: 
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and similarly 
Equation: 


N(f) =E(f) -—E(f)e?" = E(f) (1—e 8") = E(f) (ef —e Me IF = E(f)2j sin (nf)e 


and power spectrum[footnote], (use that 2 cos? (2) = 1+ cos (2a)) and 2 sin? (x) = 1— cos (2z)) 

The same formula for |M (f)| and |.N (f)|? could be obtained by computing it as the Fourier transform of the 
auto-correlation (see [link]); indeed, for m(t): r (0) = 20”, r (1) = r(—1) = o? and r(k) = 0 for all other k; 
for n(t) the same except r (1) = r(—1) = —o”. 

Equation: 


\M (f)|? = |B (f)/?4 cos? (mf) = | (f)|?2(1+ cos (2f)) = 202(1+4 cos (2nf))Rect (f) 
Equation: 
IN (A)? = |E(f)?4 sin? (wf) = |B (f)|22(1— cos (2x f)) = 202(1— cos (2nf))Rect (f) 


where o® is the total power of the original noise ¢, and S(f) = |E(f)|” = 02. Note that neither |M (f)|? nor 
|N (f)|° are flat. Verification via matlab is easy. 


Oversampling 


Oversampling 


The principle of oversampling is simple: run ADC and DAC (Analog¢-digital 
converters) at a sampling rate f, well beyond Nyquist, say at fe = Gfo where fp > 2B 
is a sufficient sampling rate. 


We call @ > 1 the oversampling rate (OSR); it is usually an integer. For audio CD, 
8 = 4 (data is sampled at f, =176'100 samples/sec but stored on the disc at fo = 
44'100 samples/sec). 


To understand the benefits of oversampling, one has to recall the typical schema [link] 
of digital signal processing (DSP). First, in the ADC and DAC parts oversampling 
reduces the losses of quality due to working too close at Nyquist rate: sampling and 
reconstruction occur with a cutoff-frequency f. far from the Nyquist rate 2B allowing 
for simple filter design or even omission of filtering. Second, as an additional benefit, 
oversampling reduces noise created by quantification, sampling or other sources of 
errors. 


Note: down- and upsampling before and after processing (DSP: e.g.: de-noising, 
detection, enhancement, storage, transmission) allows for efficient computation and/or 
transmission at low rate. 

Equation: 


band — limit ~ ADC >| M > DSP >+ M > DAC > smooth 


Noise reduction under oversampling 


In the context of oversampling, the signal has been sampled (and thereby quantized) at 
fe = Bfo, with fo > 2B and OSR 8 >> 1. Thus, we may decimate the quantized 
signal by a factor @ and still keep the full signal quality. Decimation consists of digital 
low-pass filtering at cutoff frequency 1/(28) (which corresponds to fo /2), followed by 
downsampling. 


Low-pass filtering and downsampling will not change the signal since fp > 2B, and 
thus decimation will not change the signal power (cpre. [link], [link], [link] magenta). 
However, low-pass filtering will reduce the power of the noise; the factor of reduction 
is B if an ideal filter is used. We offer several arguments for this fact. 


Spectral picture of noise reduction (ideal filter). During the first step, low-pass 
filtering, the high frequencies of the power spectrum of the noise are cut away: 


Equation: 


S(f) = P- Rect (f) > "55 ideal filter sinc + Sidea (f) = P - Rect (f@) 


After low-pass filtering, the noise is correctly band-limited (at a loss of power); 
consequently, downsampling will not change the power any further. 


Spectral picture of noise reduction (digital filter). When using a digital filter 


b = (bi,...,b,) with energy E, = b?+... +b? to filter the samples €,,, the spectrum 
will change through filtering roughly as 
Equation: 


S(f) = P- Rect (f) > np giigital filter b > Suaigita (f) ~ P- Ey - B- Rect (f£) 


Here, we used that the DFT of a well-designed digital filter b is roughly a rectangle of 
width 1/8, meaning that a fraction 1/8 of the samples bi at the center take some 
constant value A, the others are zero. It is easy to verify that A= 8 in order for 
the energy of b to be Ey, (see Comment 8). Consequently, the portion 1/8 of the k 
samples €,, will be multiplied by A, the others will be set to zero. Using [link] we see 
that S(k/K) is multiplied by A? = F,6 for —K/(28) < k < K/(28), and set to zero 
otherwise. 


Comment 8 Let us denote by q the number of samples of the DFT of b. For a well 
designed digital filter there are roughly q/@ samples equal to some constant X, the 
remaining samples are zero. Since the DFT increases energy by factor equal to the 
sample size q, we get 

Equation: 


gE, = 62+... +6 = »7q/B 


A= EP. 


According to [link] we have Ey = qP, ~ 2f. = 1/6 for a digital filter that is well 
designed, meaning that c ~ 1. Thus, such a filter achieves approximatively the same 
power reduction as the ideal filter. Typically, however, Hy is somewhat smaller than 
170. 


Noise samples 


Black: White Noise under Decimation. 
Magenta: For comparison, a signal with bandlimit B=0.0208 is subjected to the 
same manipulations as the Noise. 
Details: A = 10, P= 47/12 = 8.3: MS 5. f= 1 7 ST aa 
Ey = 0.16 (length g = 20). 


Note that the noise overwhelms the signal in rows 1 and 4; in rows 2 and 3 signal 
and noise are comparable. (The amplitude of the signal is kept unnaturally small 
for reasons of visualization). 
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Noise samples squared 


Black White Noise under Decimation. (details see figure 1) 
Magenta: Signal for comparison, 
Blue: the mean over 100 realizations of the noise. 
Red: theoretical power; dashed with an ideal filter [link], solid with a digital filter 
[link]. 


Note: the average of the squared noise terms (blue) is roughly constant at a level 
roughly equal to its power (red). 
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Estimation of S(f) by squared DFT 


Black: |DFT|?/K of White Noise under Decimation (details see figure 1). 
Magenta: Signal for comparison. 
Blue: the mean over 100 realizations of the noise. 


Red: theoretical power and S(f); dashed with an ideal filter [link], solid with a 
digital filter [link]. 
Note: the signal and its power are not changed because of band-limitation. The 
signal to noise ratio improves by a factor M when lowpass-filtering and 
downsampling (see (c)); it does not improve when only downsampling (see (d)) 
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only downsampled: 
S(f)P = 8.3 flat 


Temporal picture of noise reduction. One may argue that only the high frequencies of 
the noise are removed, which are not used anyway. In an audio signal, for instance, 
setting f, = @- 44kHz and decimating by @ removes only the noise portion in the 
frequency range which can not be heard anyway. So, no benefit should result. This 
argumentation does not take into account that the power of noise is reduced which must 
result in smaller noise values. Comparing rows 2 and 4 in [link], we note that the noise 
magnitude has indeed decreased after filtering, but not under pure downsampling. 
Comparing the spectra $(f) in [link] we note in the low frequencies from -0.1 to 0.1 
the same average level of the spectral samples of roughly P = A?/12 = 8.3; 
however, for the downsampled noise all spectral samples are at this average level, 
while the filtered noise possesses 8 = M more samples and only a portion 1/M of 
them are that level, the others are at zero. Thus, the power of the filter noise is M/Z times 
smaller than the power of the downsampled noise. (This despite the fact that only high- 
frequency noise has been removed.) 


Finally, we may also argue directly in the time domain using statistics. After filtering, 
the noise terms are nx, = 61441 + b2€K42+... bgEk+,- Intuitively, this is some sort of 
average of the original noise, since b}+..+bg = 1. But an average is always closer to 
the theoretical mean than the samples themselves, and the mean is here zero. So, 
intuitively, the filtered noise n; takes smaller values. Since the mean of the filtered 
noise terms nx is zero, a meaningful measure of their size is their power, i.e., their 
variance (cpre. [link]). 


Using simple rules from statistics[footnote] we get 
IE [n,] = bj+...+b7 IE [e] = E, - P. Short: the noise level (its power) reduces by a 
factor Ey, just as we have seen before. 


The variance of a sum of independent random variables is the sum of the variances; 
when multiplying a random variable by a constant b the variance in multiplied by b?. 


Illustration of noise reduction. 


For a illustration see [link], [link], [link]. The main comparison is row 3 (oversampled 
and decimated) with row 4 (sampled at lower frequency). While the signal is 
practically identical, the noise is reduced dramatically in row 3. 


Note also the dependence in the noise after lowpass filtering: it is possible to predict a 
few of the following noise terms. Indeed, the noise is not white and its spectrum is not 
flat. This dependence disappears almost entirely after downsampling: the prediction 
does not reach M samples into the future. 


Note the shape of the estimated S(f) is not exactly a rectangle because the filter is not 
ideal (length 25). A longer filter would result in better shapes but also in longer delay. 


Conclusion In summary, low-pass filtering reduces power by a factor 3: Only a portion 
1/ of the spectrum remains roughly unchanged the rest is set to zero. A formal 
calculation goes as follows, setting f. = 1 as usual for digital filtering and 

P = A?/12 for a quantization noise (why we are interested in noise, and not in the 
noisy signal, see Comment 9) 

Equation: 


1/2 g 1/(28) 1 1 A? 
Pea iltere df = Pdf = —P = —— 
20 a —1/(28) : B B 12 


Comment 9 We do not compare power of signal to power of noisy signal, because they 
are almost equal; also, a theory of the power of the noisy signal is hard to develop since 
power(signal+noise) is not power(signal)+power(noise), unless signal and noise are 
uncorrelated, e.g. if the signal changes slowly and the noise has mean 0. In conclusion, 
we compare the power of the signal to the power of the noise. The ratio is called SNR 
“Signal to Noise Ratio”. 


Since downsampling of a properly band-limited signal does not change its power (see 
[link], [link], [link]), the factor 1/8 is the overall change of power during decimation. 
In other words, the SNR[footnote] improves under oversampling with subsequent 
decimation by 8 (which includes low-pass filtering at 1/28!) by the factor £, i.e., it 
improves by the oversampling rate (OSR). For a statistical argumentation, see 

box Comment 10. 

The Signal to Noise Ratio (SNR) is the ratio of the power of signal and power of noise. 
Using Parseval's equation, it can be computed in the time or frequency domain: 


Equation: 


x(t) “dt  ——-X(f) “df 


SNR = 
n(t) “dt = -N(f) “df 


Comment 10 Statistical explanation of noise reduction. Short: filtering is a special way 
of averaging the samples. Averaging reduces variance, thus it reduces the noise power. 
In addition, the averaging is done in a special way as to leave small frequencies, i.e. 
slowly varying components, intact. Thus, the averaging reduces only noise, which is 
high frequency, i.e. quickly varying, and leave the signal intact. Somewhat more 
rigorously: Filtering noise means to compute 

Equation: 


M 
/ 
é, = bnEk—n 
n=1 


Computing the variance of the resulting, filtered noise, we find: 
Equation: 


M M 
a 2 — 2 
var €, = Ob var(ez,_n) =o be. 


Here, we use that the original noise is independent. See footnote [link]. Now, for a FIR 
filter b it is easy to verify, e.g. with matlab, that ae b2 ~ 1/8. This approximation 
becomes more accurate the larger / is; it reflects the fact that the power of bis close 
to the power of the ideal lowpass filter c and the power of c is exactly 1/8. A similar 
approximation holds also for any least square filter a: This is due to the fact that a 
approximates the ideal lowpass filter c in the least square sense. More precisely, let us 
denote the norm of ¢, i.e. the “length” of c by || ¢ ||. A well known formula says that 


| ¢ |= cf +... +c}. To approximate c in the least square sense means that 


|| a — c || is made as small as possible. But another well known fact says that 

| || a || — || e || | <|| @—c |]. This means, that the norms of a and c, thus their 
energies, are almost identical. In summary, 

Equation: 


var e, —o°*/8. 


In other words, the average power per sample is reduced be a factor (. 


Noise-Shaping 


Noise-shaping 
Noise shaping is a procedure which allows to put some color on an additive white 
noise. While this adds to the overall noise, one can shape the resulting noise such that 


its prevalent “color” lies in the high frequencies and at the same time reduce the 
presence of noise in the low frequencies. 


E,(f) Xo(f) 


Usual quantization noise. 


When Oversampling, only the part in the band [— fy /2, fo /2] will be 
relevant after down-sampling. 


Noise Shaping moves the noise away from the smaller frequencies. 


In the context of oversampling, noise shaping achieves an improvement if it is done so 
that it reduces noise in the band [— fp /2, fp/2] and pushes it to the “colors” with 


frequencies in |fo/2, f./2] and [—f./2, —fo/2]. 


The most simple version of noise shaping results in two noise terms, the one produced 
by the system, one added with a delay on purpose. Set €g = O, then: 
Equation: 


Lk —> —EK-1 —> (ay = Ek-1) => 


Q i Uk output 
> =" 
oe Ek = Yk — (fe —Ex-1) feedback 


This assumes that we know or can know the error. Such is the case with quantization, 
where we can compute €, simply as the difference between output and input of the 
quantization (see [link]). This is not feasible in the example of a wireless channel with 
inference. 


To study the effect of noise shaping we compute the Fourier transform of the overall 
eIror ny = YR — LE. We find 
Equation: 


N= Yk — Lk = Ek — Ek-1 


or 
Equation: 


with Fourier transform (use that f. = 1/T~) 
Equation: 


E(f)—E(f)e#""f = E(f) 1-7? lfe 
E(f)2 sin (wf/fe)e "1!" 


| 


N(f) 


| 


and power spectrum, (use that 2 sin? (2) = 1— cos (22)) 


Equation: 
IN (A)? = 1B (/2(1- 00s (24) 


e 


Note that the spectrum is no longer flat; in other words, the noise nx is colored. We 
note that the colored spectrum is small for small f, and it is still spread over a period 
of length f-. Thus, noise shaping results in a reduction of noise in the small 
frequencies (see [link]). 


Now we continue as with simple oversampling: since the signal has been sampled at 
fe = Bfo, with fo > 2B and GB > 1 we may low-pass filter after noise shaping at 
cutoff frequency fo /2. Thus, we take advantage of the fact that most of the power of 
the “shaped” or “colored” noise is in the “‘high'-frequency bands [fo /2, f-/2] and 
[—fe/2, —fo/2] with very little noise at frequencies around 0. 


To assess the gain we compute the power of the noise after this low-pass filter. Note 
that the low-pass filter sets N(f) = 0 for f./2 > |f| > f/2: this removes much of 
the noise as demonstrated in the following computation; it also guarantees that power 


won't be changed when downsampling after filtering to fp. Using again that 
E(f) = 42/12 we get: 
Equation: 


1 fo/2 


P. shaping — 


fo/2 anf 
ia N(f)?df = —2 7 *2(1- ( )) 
fe se i) j fe 0 | Me “ I 
_ 4a fe .. ( 2mf\\ fo/a _ a) : ~ sin (7 
= ag (f~ 8 (F)) 8 = Sel ae - ae a (F)) 


Using the approximation sin (x) ~ x — 23/6 + 2°/5! — +... we find 
Equation: 


Pshaping =) 1 1, Tv oe 1 _ 1 Tv _ 13 7 12 “8 
Paw “(a5 ~ ae (5) =*(Ca5- ae (F- e)) =F 


Noise reduction under oversampling with noise-shaping 


In conclusion, the SNR improves (as compared to no oversampling) under noise 
shaping by the inverse of the factor [link], thus roughly by + G3. In decibel, this 
corresponds to approximatively 10 log) (3) — 20 logy, (7) + 30 log,, (8), thus, 
roughly 3 times more than with oversampling alone. 


Numerical Examples: 
With an OSR of @ = 128 the SNR improves 


¢ under oversampling alone by 10 log;, (128) ~ 21dB 

e under oversampling coupled with basic noise-shaping as above by 
approximatively 10 log,,) (3) — 20 logy, (7) + 30 log,, (128) = 58. 0445dB, 
or roughly 3 times more. 


When doubling the OSR the SNR improves 


¢ under oversampling alone by 10 log), (2) ~ 3.01dB 
¢ under oversampling with noise-shaping by approximatively 30 log,,) (2) = 9.03 
dB, or approximatively 3 times more. 


